Mechanically-driven growth and competition in a Voronoi model of tissues

The mechanisms leading cells to acquire a fitness advantage and establish themselves in a population are paramount to understanding the development and growth of cancer. Although there are many works that study separately either the evolutionary dynamics or the mechanics of cancer, little has been done to couple evolutionary dynamics to mechanics. To address this question, we study a confluent model of tissue using a Self-Propelled Voronoi (SPV) model with stochastic growth rates that depend on the mechanical variables of the system. The SPV model is an out-of-equilibrium model of tissue derived from an energy functional that has a jamming/unjamming transition between solid-like and liquid-like states. By considering several scenarios of mutants invading a resident population in both phases, we determine the range of parameters that confer a fitness advantage and show that the preferred area and perimeter are the most relevant ones. We find that the liquid-like state is more resistant to invasion and show that the outcome of the competition can be determined from the simulation of a non-growing mixture. Moreover, a mean-field approximation can accurately predict the fate of a mutation affecting mechanical properties of a cell. Our results can be used to infer evolutionary dynamics from tissue images, understand cancer-suppressing effects of tissue mechanics, and even search for mechanics-based therapies.


I. INTRODUCTION
Tissues are communities of cells and undergo evolutionary processes.Cell competition, in which one cell type takes over another, is particularly relevant for cancer invasion, but also in developmental processes [1][2][3].Cell competition between a mutant and a resident population is the result of changes that alter the growth rate of the mutants with respect to the resident.If a mutant acquires a higher growth rate, it will eventually fixate and take over the entire population.
In normal tissues, there is a turnover of cells over time [4]: cells undergo division and are eliminated by apoptosis or cell extrusion [5].The balance between cell division and extrusion controls the overall growth of the tissue.In what follows, we will refer to cell extrusion without distinction between apoptosis and extrusion because, at our level of description, they both contribute to the removal of a cell from the monolayer.We are interested in this paper in the mechanical regulation of growth and its implications in cell competition [1,6,7].At the cellular level, there are mechanotransductive proteins that provide mechanical feedback to cells, whose levels can promote or hinder growth [8][9][10].On a macroscopic level, it has been shown that exerting pressure on cell spheroids hinders growth, both by lowering the division rate and by increasing the rate of apoptosis [11].In cell monolayers, cells try to maintain a constant homeostatic density by extruding cells in denser areas, where the pressure is greater [12].The idea that cell growth is regulated by pressure is not new [13,14], but to our knowledge its implication for cell competition has not been explored outside of continuum models [15].
To investigate the effect of mechanics on cell competition, we study a commonly used model of a confluent cell monolayer: the Self-Propelled-Voronoi (SPV) model.We supplement it with a stochastic growth process with birth and death rates controlled by tissue mechanics.In the SPV model, the motion of the cell centers is derived from an energy functional, and the network of cells is constructed using a Voronoi tessellation.The energy functional is the same as in the Vertex Model [16,17], but we model the motion of the cell centers rather than the vertices.An important aspect of both models is that they predict a jamming-unjamming transition [18][19][20] that has been suggested to play a central role in several biological processes [21][22][23].
After introducing the SPV model for a heterogeneous population composed of two types of cells and its birth-death dynamics, we simulate invasion scenarios with all possible singleparameter mutants to determine the outcome of the competition.We find that the fate of a mutant can be predicted from a simulation of a mutant-ancestor mixture without cell division or death.Using this fact and a mean-field approximation, we explain the results of our simulation by a simple argument based on energy minimization.

A. Heterogeneous Self-Propelled Voronoi model
We use the Self-Propelled Voronoi (SPV) model to describe the dynamics of a confluent monolayer of active cells, as found in tissues [19].In this model, the cell centers evolve A time unit represents 100 simulation time steps.according to equations of motions derived from an energy functional, with self-propulsion motility.Once the positions of the cell centers are computed, the shapes of the cells are determined using a Voronoi tessellation, i.e. by maximizing the size of each cell with respect to the position of its cell center.Other popular models such as the Vertex Model [16] or the Cellular Potts model [24] use more degrees of freedom to compute the shape of cells, but the simple aspect of having the cells themselves as degrees of freedom make evolutionary dynamics much easier to track, especially since we have little interest in the details of the cell shapes.
The energy functional for the cells is the same regardless of the chosen model: where the index α represents a cell in the tissue.Although the requirement for the cells to form a Voronoi tessellation appears as an added constraint compared to the Vertex Model, the dynamics of the network of cells is very similar in both cases far from the low noise limit [25,26].
We are interested in the competition between two types of cells, a resident population labeled with a subscript r and a mutant population labeled with a subscript m .We assume for now that cells are homogeneous within each type.We define the total number of cells N = N r + N m , and the total energy of the population of cells is: A common choice is to normalize all lengths by A 0 r , fixing the normalized preferred area to unity.However, because our simulation using cellGPU [27] is setup to have a fixed volume L 2 set by the initial number of cells, we believe it is appropriate to normalize the lengths by L/ √ N 0 , where L is the length of the simulation box and N 0 is the initial number of cells.
With this normalization, the average cell area is determined by the cell number at time t: The energy is normalized by K r L 4 /N 2 0 .The dimensionless total energy e reads There are two dimensionless variables, the areas a α = A α N 0 /L 2 and the perimeters p α = P α √ N 0 /L, and 7 dimensionless parameters.Among those 7 parameters, four are so-called geometrical parameters a 0 r , a 0 m , p 0 r , p 0 m and three are mechanical parameters K, Γr , Γm .At each timestep, the force exerted on cell center i is given by where NN(i) refers to the nearest neighbors of cell i.The equation of motion for cell center i is with µ the mobility of cells taken to be homogeneous amongst cells, and n i = (cos θ i , sin θ i ) a unit polarization vector that randomly changes direction over time: The cell centers are displaced using a forward Euler discretization of eq. ( 6) and the Voronoi tessellation is computed with the new position of the cell centers.

B. Coupling between mechanics and growth
While we can infer general trends for the effect of mechanics on growth such as pressureinduced inhibition of growth [11] or regulation of homeostatic density [12], pinpointing the exact influence of mechanical variables on individual cells proves challenging.Experimentally, researchers have access to macroscopic measures such as the force on a cell assembly or the force exerted on the substrate using traction force microscopy.However, quantifying bulk properties such as the stress within a cell monolayer remains elusive.The stress within a monolayer can be computed using imaging-based inference methods [28,29], but since the methods rely on formulating models similar to the one we are using, it does not provide new information.
In the absence of a clear coupling between mechanics and growth dictated by experiments, we can think of a generalized linear response of growth to mechanical changes.If one thinks of the mechanical parameters a 0 , p 0 in the energy functional as some kind of target homeostatic value for the tissue, then the birth and death rates of a cell r i can be expressed as The first term corresponds to the basal components of the birth and death rates, respectively r b i 0 and r d i 0 , which do not depend on mechanical parameters.The second term is a linear coupling between the mechanical state of the system and the growth rate r b i − r d i , with a coupling coefficient λ and a parameter β that weighs the contribution of the perimeter with respect to the area.The min and max in the death and birth rates ensure that each rate is positive.Higher coupling terms are neglected.
In this first work, we focus on the case β = 0, which means that the growth rate depends only on the area of the cells, or equivalently on the density in the tissue.This particular coupling is relevant in epithelial monolayers, where there is more extrusion of cells in denser regions [12].This coupling is related to the hydrostatic pressure defined as The hydrostatic pressure represents the normal forces exerted on the edges of the cells [29].
For simplicity, let us consider the case with the same basal birth and death rates for all cells, r b i 0 = r d i 0 = r 0 .In that case, the rates are given by: With our normalization scheme, the average growth rate at time t is fixed by eq. ( 3) and depends on the number of cells: The steady state is reached when the average growth rate ⟨g i ⟩ vanishes and the steady-state cell number is set by the initial cell number and the preferred area N SS = N (t = 0)/a 0 .We represent the evolution of cell number over time using different growth couplings in fig. 2.
Another mechanical observable associated with growth is the pressure [11,14,30], which corresponds to a special value of β.The pressure is defined as the opposite of the isotropic part of the stress tensor, which can be determined from the forces using the virtual work principle [31].The dimensionless pressure P, normalized by KL 2 /N 0 , is given by: To linear order in a − a 0 and p − p 0 , it corresponds to β = Γp 0 /(2a 0 ).The birth and death rates of cell i are expressed as: where P h is viewed as the homeostatic pressure, since it determines the value of the pressure in steady state.To keep the number of initial cells constant, there is a unique value of P h that depends on the mechanical parameters.
In a future work, a comparison between different couplings between the growth rate and mechanical properties and their effect on cell competition will be studied.Notably, the limit β → ∞ is unstable.When β ≫ 1, a larger p leads to a higher growth rate, regardless of the area.Since there is no constraint on the perimeter, this leads to infinitely elongated cells.

A. Growth of an homogeneous tissue
The cell number over time from different couplings between mechanics and birth-death rates is shown in fig. 2. In all cases, r b i 0 = r d i 0 for all cells.We looked at three cases, λ = 0 (no coupling to mechanics), β = 0 which is density regulation, and coupling to pressure.In each case, there are five realizations of the stochastic birth-death process using a Gillespie algorithm.We see that when there is no coupling to mechanics, the stochastic birth-death process does not reach a steady state.The cell number is only subject to demographic noise and there is no effective carrying capacity.However, coupling to mechanical forces stabilizes the number of cells to a value that depends on the mechanical parameters.We therefore argue that a mechanical feedback is necessary to maintain a fixed average cell number when the birth-death process is noisy.

B. Invasion by a mutant population
After rescaling physical units, our model contains a total of 15 parameters: 7 energy parameters, 4 growth rate parameters, and 4 motility parameters.Since we are interested in the competition between mutants and residents, we decided to fix some parameters of the resident population.We end up with a total of 8 parameters involved in the competition.
The parameters and the effect of their change are recapitulated in table I.
The parameters of resident and mutant cells can be thought of as phenotypic traits that can be changed by mutations.We are interested in the evolutionary outcome of a change in any of these phenotypic traits.To test this, we performed numerical invasion experiments in We notice that in the solid-like state (solid line) there is a negative pre-stress value that is not relaxed.However, the liquid-like state (dashed line) is able to relax this pre-stress.
which we introduce a small fraction f of mutants among a resident population.The mutants have a small change of the order of 10% in one of the parameters compared to the residents, and we performed 8 invasion experiments, one for each of the competition parameters (see table I).Looking at the evolution of the fraction f of mutants over time, we can determine if the change is advantageous -f increases on average, eventually leading to fixation -or deleterious -f decreases on average leading to extinction.The evolution of the ratio between the mutant fraction and the resident fraction is represented in fig. 3 on a semi-log scale.On average, the mutant fraction evolves according to a logistic equation

TABLE I: Simulation parameters and their effect on fitness
The effect of the basal rate of the mutant r 0 m is trivial since it provides a constant growth rate difference.For the two geometric parameters, the preferred area a 0 and the preferred perimeter p 0 , it could be expected that the growth rate difference generated would change as the mutant fraction evolves.However, we find that the growth rates of the resident and mutant population change over time, but the difference between the two remain constant, mainly due to the volume constraint.Let us take for example the case where the preferred area is varied, and the mutant has a smaller preferred area a 0 m = 0.9 while the resident preferred area is a 0 r = 1.The volume constraint fixes the average area of the entire population by the number of cells, according to eq. (3).Initially, the average area of each cell is fixed to unity.The mutant cells introduced are trying to be smaller than the resident cells to minimize their energy.With the average area fixed roughly to unity, the mutant cells will be greater than their preferred area, but so will the resident cells that grow in size to compensate for the smaller mutant cells.This effect increases the growth rate of both mutant and resident cells, depending on their relative fraction and the total number of cells.
However, the difference between the two growth rates does not depend on the fraction nor the total number of cells.A more detailed derivation from energetic arguments is provided in section III C. The argument for a change in the preferred perimeter p 0 follows the same path, since it effectively changes the realized area of the cells.We find that the effect of the geometric parameters a 0 and p 0 on fitness is stronger in the solid-like phase than it is in the liquid-like phase.We attribute this effect to the broader distribution of cell areas in the solid-like phase with respect to the liquid-like phase (see fig. 4).
A small change in the strength of the coupling between the area and the growth rate λ seems to have no effect on the evolution of the population.Mechanical elasticities Γ, K play little role in the fitness when they are varied alone.
However, we show in section III C that they control the magnitude of growth rate difference, even though they cannot create such a difference.

C. Evolutionary outcome determined by a non-growing mixture
Since the rate of change of the mutant fraction is constant over time, we can determine the fate of the competition between mutant and resident strains by simulating a mixture of the two strains at a constant fraction, without growth.We simulate a 50-50 mixture of mutant and resident strains without growth and compute the average realized areas for mutant and resident strains, resp.⟨a m ⟩ and ⟨a r ⟩, which gives the average growth rate difference between mutant and resident strains s = λ(⟨a m ⟩ − a 0 m − (⟨a r ⟩ − a 0 r )).The results are shown in fig. 5.The sign of the growth rate difference is well predicted by simulating a mixture without growth, but the actual value of the growth rate is somewhat different.This discrepancy could be explained by the fluctuations in cell number induced by the fluctuations in the birth-death process, since a change in cell number changes the average value of the area (see eq. ( 3)).Experimentally, this mean the evolutionary information can be obtained at any point in time from an heterogeneous population of cells, by averaging the distribution of areas.
In this section, we present an energy argument to explain why the growth rate difference between mutants and residents does not depend on the fraction.Consider a mixture with residents having a preferred area a 0 r = 1 and mutants having a preferred area a 0 m .Let us assume that both residents and mutants have the same preferred perimeter p 0 .The initial when the birth and death rates are controlled by the area, as stated in eqs.(10) and (11).
Each panel corresponds to an average of 10 identical simulations, where there is a small change in one parameter for the mutant population.Solid lines are for simulations where the resident is in a solid-like state with p 0 r = 3.3, and dotted lines are for simulations where the resident is in a liquid-like state with p 0 r = 4 population is composed of N 0 cells in a fixed volume L 2 (an average cell initially has an area of value 1).We consider that at time t there are N t cells, a fraction f = N m /N t of mutant cells, and 1 − f = N r /N t resident cells.Consider a mean-field approximation to solve the evolution of identical average cells.In that case, the volume conservation at time t reads: FIG. 4: Distribution of areas and perimeter for solid (upper row) and liquid (lower row) 50-50 mixtures of resident and mutants averaged over time.In the left column, mutants have a preferred area of a 0 m = 0.9 while the one of residents is a 0 r = 1 (horizontal dashed lines).All other parameters are equal.In the right column, mutants have a preferred perimeter of p 0 m = 0.9 while residents have a preferred perimeter p 0 r = 1 (vertical dashed lines).All other parameters are equal and have their reference value.Equation ( 16) sets the average cell size to be N 0 /N t .However, this value does not minimize the energy of the mixture with two different preferred areas.We assume that there are small deviations δa m , δa r ≪ N 0 Nt w.r.t the average value, and write the areas as a m = N 0 /N t + δa m and a r = N 0 /N t + δa r .If the deviations in areas are small, to first order the perimeters are We can make the assumption that the proportionality coefficient c is the same for both mutants and residents since having the same preferred perimeter means they are similar shaped polygons.mutant is a 0 m = 0.9 while the one of residents is a 0 r = 1.The solid line represents the simulation of the birth-and-death process in a solid state p 0 r = p 0 m = 3.63.All other parameters are equal and have their reference value.The dotted line is the estimation from the mean-field approximation given by eq. ( 19) and the dashed line is an estimation using the area difference given by a mixture without growth.We see that although both estimation predict the same evolutionary outcome, they overestimate the the rate of evolution.With this perturbation, there is only one degree of freedom, and we can choose it to be δa m .By minimizing the energy with respect to δa m , we obtain the mean-field equilibrium areas for the mutant and resident cells.From these areas, we determine the mean-field prediction for the fitness advantage when the preferred area is varied: with ∆a 0 = a 0 m − a 0 r is the difference in preferred areas.The sign of ∆a 0 determines the evolutionary outcome and supports the finding from the invasion simulation that a 0 m < a 0 r leads to a fitness advantage.
The exact same thing can be done when the preferred perimeter is varied, and gives with ∆p 0 = p 0 m − p 0 r is the difference in preferred perimeters.Similarly, the sign of ∆p 0 determines the evolutionary outcome and supports the finding from the invasion simulation that p 0 m > p 0 r leads to a fitness advantage.Both eqs.(19) and (20) confirm that the fitness is independent of frequency.The growth rate of each type depends on the frequency of the mutant, but the difference between the growth rate of each does not depend on the frequency.
We can also note that although the normalized perimeter elasticity Γ does not play a role without a mutation in the preferred area or perimeter, it controls the rate of evolution when such a mutation exists.The greater Γ is, the faster evolution occurs.A hand-wavy argument is that a high perimeter elasticity fixes the shape of the cells, making it harder to deform in the network to accommodate changes in area.
Both rates depend on the coefficient c that relates a change in area to a change in perimeter.For a regular n-sided polygon, there is a relationship a = 1/c n p 2 , with c n = 4n tan(π/n).In the solid-like phase, when averaging over time the distribution of areas and perimeters, we indeed observe a linear relationship between the area and the squared perimeter, as seen in the first row of fig. 4. Using the coefficient for regular hexagons c 6 ≈ 13.9 fits well the distribution of areas and perimeters in the solid phase.In the liquidlike phase, the shapes of the cells are not well described by regular polygons, and there is no a priori relationship between the area and the perimeter.The simulation of a mixture of resident and mutant cells shows that there is an area change induced by a change in preferred perimeter and vice versa.This change in area is smaller than for the solid-like phase, which explains the smaller rate of change of the mutant fraction in the liquid-like phase compared to the solid-like phase observed in fig. 3.

IV. DISCUSSION
Mechanics play a role in the regulation of growth in tissue [11,12,30], and thus can be a driver of cell competition, for example in the case of cancer invasion [7].Using an agent-based model of tissue mechanics with mechanically-regulated growth bridges the gap between evolutionary models of tumors [32,33] that are concerned with the genetic diversity in tumors and macroscopic hydrodynamic or elastic models that look at deformations and growth induced by heterogeneous cells [34,35].
Mechanically-regulated growth in a homogeneous tissue controls the average number of cells in the system even in the presence of demographic noise, by effectively providing a carrying capacity.To study competition in heterogeneous tissues, we simulated invasion scenarios by different mechanical mutants and have identified that the geometrical parameters in the Voronoi model, the preferred area a 0 and the preferred perimeter p 0 , confer a fitness advantage to the mutants, while other parameters do not seem to play a role.
A change in preferred area can be interpreted macroscopically as a change of homeostatic pressure, where lowering the preferred area is similar to increasing the homeostatic pressure, i.e. the pressure at which the tissue is in a homeostatic growth state.In this light, this model agrees with hydrodynamic models of tumor growth coupled to mechanics [13,14].
The fitness advantage, when it exists, is constant over time, which means that it does not depend on the frequency of the mutant in the tissue.Therefore, we can explain the long-term evolutionary fate of a mutant -fixation or extinction -either by energy consideration on a tissue made of average cells or by simulating a mixture of mutants and residents without growth.This is relevant for experiments since it means that one does not have to conduct a full evolutionary experiment to infer the outcome, it is sufficient to know the average state of tissue at a certain time point.
In this work, we studied in detail the case in which the growth rate depends on the cell area, and we should note that the results obtained depend on the chosen mechanical coupling.For example, if one were to consider the birth rate of a cell to be proportional only to the realized area of the cell instead of the difference between the realized and preferred areas, a change in preferred area would have opposite effects on the mutant fraction, as a larger preferred area would just lead to larger cells.A systematic study in future works of possible couplings between mechanics and growth could help discriminate which coupling is most relevant to experimental realizations.We found that one of the key features of the SPV model, the jamming-unjamming transition [18][19][20], is not a key driver of cell competition.
We believe that there are two reasons why there are no fundamental differences between the two phases.The first is that the growth process essentially fluidizes the tissue [36] and gives rise to a non-vanishing effective diffusion, effectively unjamming the jammed phase.
The second reason is that over the time of the simulation the system is essential well-mixed, there is no sustained spatial difference between the mutants and the residents.A next step for this work would be to study the effect of mechanics on genetic diversity.
At larger spatial scales, one would observe nontrivial growth-induced pressure gradients and traveling fronts, which could change important aspects of the competition [15,37,38].

FIG. 1 :
FIG. 1: Overview of model.a) A Voronoi cell and its parameters.b) Example of a cell division.The red cell is split in half by a segment that goes through the center of the cell, with the angle of the segment chosen at random.Two new cell centers are created at a certain distance normal to this segment, and a new Voronoi tessellation is computed.c) Example of cell death/extrusion in which a cell center is removed and a new Voronoi tessellation is created.d) Initial and e) final configuration of an invasion simulation, where a small droplet of mutant with a lower preferred area is introduced with fraction f = 0.1.

FIG. 2 :
FIG.2: Left panel: evolution of the number of cells over time for different regulation of the birth-death process.We notice that purely stochastic process (blue) leads to large deviations of the cell number over time, while birth-death processes regulated by mechanics, either through the density (red) or the pressure (green), lead to a steady-state value of the number of cells.Right panel: distribution of the time-averaged pressure in the tissue.We notice that in the solid-like state (solid line) there is a negative pre-stress value where s = g m −g r is the growth rate difference between mutants and residents.If s does not depend on the mutant fraction, the solution to the logistic equation is simply ln(f /(1 − f )) = st.Therefore, the straight lines on figure 3 correspond to growth rate differences that do not depend on the mutant fraction.From figure3we identify 3 parameters that provide a fitness difference, while the other 5 parameters do not significantly change the fraction of mutant.Interestingly, for the 3 parameters that provide a change in fitness, this change is constant over time, since ln(f /(1 − f )) is a linear function of time.Strength of mechanical coupling to growth No apparent effect.Resident value fixed to λ r = 1 r 0 r,m Intrinsic growth rate Provides a fitness advantage when r 0 m > r 0 r .Resident value fixed to r 0 r Resident value fixed to D r = 1

FIG. 3 :
FIG.3: Evolution of the fraction of mutants from an initial droplet of fraction f = 0.1,

FIG. 5 :
FIG. 5: Evolution of the fraction of mutants over time when the preferred area of the